Gnuplot est un logiciel pour le tracé graphique de fonctions (2D et 3D) et l'analyse de données. Certes, gnuplot est moins "glamour" que Desmos, mais il est purement local, et pour traiter des fichiers de données, il est le meilleur.

Mais gnuplot est aussi un calculateur interactif en mode ligne de commande (terminal). Il suffit de taper "print" suivi d'une expression, ou "plot" suivi d'une fonction. La liste de fonctions mathématiques prédéfinies est impressionnante, (elle contient autre autres les fonctions de Bessel...) mais on peut en définir et ajouter d'autres !

J'ai ajouté certaines fonctions qui me paraissaient utiles, ou simplement très astucieuses, dans le fichier GNUPLOT.INI qui est sur mon ordi dans D:\Documents\installateurs de programmes :

Notons que toutes ces définitions de fonction listées ci-dessous utilisent la récursion de préférences aux boucles. Ce n'est que très tardivement que j'ai découvert que gnuplot acceptait les syntaxes do for [i = 1:5] {pr i} et i=10; while(i>0){i=i-1; pr i}dans les commandes mais malheureusement pas dans la définition de fonctions. Par contre gnuplot accepte des blocks de fonctions : taper help function blocks
exemple

On aurait pu écrire plus simplement $sinc(x)=(x==0)?1.0:sin(x)/x, mais c'est un exemple.

J'ai noté aussi une fonction utile pour créer des arrays de taille dynamique : split("a b c") renvoie ["a","b","c"]. Pour la taille des arrays on utilise l'opérateur infixé |...| Par exemple |split("a b c")| renvoie 3

donne le plus petit facteur premier qui divise floor(x)(1 si premier)
facto(n) =facto2(floor(n))
facto2(n) =(n<=1)?1:(n%2==0&&n!=2)?2:(n%3==0&&n!=3)?3:(n%5==0&&n!=5)?5:(n<49)?1:facto1(7,n)
facto1(d,n)=(n%d==0)?d:(d*d>n)?1:facto1(nextprime1(d),n)

donne le prochain nombre premier après floor(x)
nextprime(n) =nextprime2(floor(n))
nextprime1(n)=(facto2(n+2)==1)?n+2:nextprime1(n+2)
nextprime2(n)=(n<2)?2:(n%2==0)?nextprime1(n-1):nextprime1(n)

donne tous les facteurs sous forme de chaine eg 24->" 2 2 2 3"
factorise(n)=(facto(n)==1)?" ".n:" ".facto(n).factorise(n/facto(n))

Noter que split(factorise(12)) renvoie l'array ["2","2","3"] ; c'est un array de strings, mais split(factorise(12))[3]+2 renvoie quand même 5 (dans gnuplot le premier index des array est 1, pas 0 comme en js ou Python)

Donne le nième nombre premier, en commençant par nthprime(1)=2, nthprime(2)=3etnthprime(0)=1, ce qui est assez logique
voir aussi https://t5k.org/nthprime/index.php#nth

Version très limitée ; (n<709) à cause des limites de récursion de gnuplot :
nthprime(n)=(n<=0)?1:(n==248)?1571:(n==479)?1407:nextprime(nthprime(n-1))

Version sans limite ?
Attention, au delà de n=10000, le temps de calcul peut devenir assez long (4 secondes pour 15000)
nthprime(249)=1559 ; nthprime(10000)=104729
nthprime(n)=(n==0)?1:nthp(0,n)
Explication :
La fonction suivante donne le kième nombre premier après n (si k=0, renvoie n)
nthp(n,k)=(k<=0)?n:(k==1)?nextprime(n):nthp(nthp(n,k/2+k%2),k/2)

pgcd de deux entiers (algo d'Euclide) : Plus Grand Commun Diviseur
pgcd(a,b)=pgcd1(floor(a),floor(b))
pgcd1(a,b)=(a<b)?pgcd(b,a):(b==0)?a:pgcd1(b,a%b)

Plus Petit Commun Multiple
ppcm(a,b) = a*b/pgcd(a,b)

multiplication binaire sans retenue. cf fichier factomul.txt
par exemple mul(5,5)=17
mul(a,b) =mul1(floor(a),floor(b))
mul1(a,b)= (b<2)?(b&1)*a:((b&1)*a)^(mul1(a,b/2)*2)

Facultatif : R(a,b)=floor(a)*floor(b)-mul(a,b) ->retenues de la multiplication

g(n)=floor(n)^(floor(n)/2)

L'intégrale de parité, chère à Gérard Langlet, telle que g(ig(x)) = ig(g(x) = x où g(n) est le code gray
ig(x) = ig1(floor(x))
ig1(n) = (n<=1)?n:(ig1(n/2)%2)^(n%2)^(ig1(n/2)*2)
une version plus rapide : le nombre de récursions est le nombre de "1" dans l'écriture binaire de x. C'est encore plus rapide que la FFT !
ig(x) = ig2(0,floor(x))
ig2(r,n) = (n<=0) ? r: ig2(r^n, n/2)

en commençant par F(0)=0, F(1)=1 donc F(2)=1, F(3)=2, F(4)=3, F(5)=5 etc.
c'est une version modifiée de Fast doubling Fibonacci algorithm (JavaScript) (page web) (bien plus rapide que l'approche naïve)
Pour les négatif, on utilise les negafibonacci : F(-4) = -3, F(-4)= +5. Cf https://en.wikipedia.org/wiki/Zeckendorf%27s_theorem
On utilise la représentation complexe : fib(n) renvoie le complexe {F(n), F(n+1)}

Notons qu'on aurait pu utiliser la formule de Binet avec , mais la fonction gnuplot ci-dessus a l'avantage de ne manipuler que des entiers (de Gauss, certes...), et surtout de fonctionner pour les négatifs.

Noter, ce qui est peu connu, que et que

voir aussi représentation de Zeckendorf de n

conversion de n decimal vers base b<= 10 :
d2b(n,b)=(n<=1)?n:n%b+10*d2b(n/b,b)

conversion base b (<=10) vers décimal
b2d(n,b)=(n<=1)?n:n%10+b*b2d(n/10,b)

conversion de n de la base d (départ) vers la base a (arrivée). Les deux bases doivent être <=10
cbase(n,d,a)=d2b(b2d(n,d),a)

Z(n) donne la décomposition "canonique" en somme de nombres de Fibonacci de n
il faut avoir implémenté F(n) auparavant, c.f. suite de Fibonacci
cf Efficient algorithm for Multiplication of Numbers in Zeckendorf Representation (page web)
e.g. Z(6) = 100100_2 = 20, bits(Z(30)) = 101000100

rend sous forme de chaîne les bits du nombre n
exemple bits(11) = '1011 '. Rappelons que dans gnuplot le . est la concaténation, comme en php
bits1(n) = (n==0)?"0":(n==1)?"1":bits1(n/2).bits1(n&1)
bits(n)=bits1(floor(n))." "
voir http://gnuplot.sourceforge.net/demo_5.2/stringvar.html pour des exemples de manipulation des chaines

rend les n premiers bits du reel x en base 2, avec le point à la bonne place
e.g. nbits(pi, 10) = 11.00100100
nbits(x,n) = bits1(floor(x)).".".nbits1(x-floor(x),n-nbb(floor(x)))
pour ci-dessous, x doit être entre 0 et 1
nbits1(x,n) = (n<=0)?"":(floor(x*2)==0)?"0".nbits1(2*x,n-1):"1".nbits1(2*x-1,n-1)

crée un nombre dont les bits pairs sont ceux de b et les bits impairs sont ceux de a. Donc un nombre formé à partir de deux autres dont les bits sont entremêlés. Cela a un lien avec les courbes de Lebesgue et de Morton.
cf 'facto par insertion de zeros.txt'
ins(a,b) =ins1(floor(a),floor(b))
ins1(a,b)= (a==0&&b==0)?0:(b&1)+2*(a&1)+4*ins1(a/2,b/2)

donnent le nombre formé uniquement des bits pairs, resp impairs
pour tout n, ins(impairs(n),pairs(n)) == n
pairs(123) = 13
impairs(a) =impairs1(floor(a))
impairs1(a)= (a<=1)?0:((a/2)%2)+impairs1(a/4)*2
pairs(b) =pairs1(floor(b))
pairs1(b)= (b<=0)?0:(b%2)+pairs(b/4)*2

Rend le b-ième bit de n (n entier, bit 0=LSB)
nthbit(11,0)=1, nthbit(11,1)=1, nthbit(11,2)=0, nthbit(11,3)=1
nthbit(n,b)=((n&(1<<b))>0

nbb(6)=3
Pas forcément le plus rapide.
nbb(x) = floor(log2(x))+1 # ok si x >= 1
version alternative plus rapide (et ok pour 0 : nbb(0)=1)
nbb(n) = (n<=1)?1:1+nbb(n/2)

nombre de bits "1" dans l'écriture binaire de n.
eg nbb1(14) = 3, nbb1(15) = 4, nbb1(0) = 0
nbb1(n) = (n<=0)?0: 1+nbb1(n & (n-1))

Inverse les bits i et j de n. (bit 0=LSB).
bitswap(11,2,1) = 13 parce que 1011 -> 1101
A noter que cela peut éventuellement rendre un nombre plus grand que n si i ou j sont plus grand que le nombre de bits de n
bitswap(n,i,j)=(((((n>>i)^(n>>j))&1)<<i)|((((n>>i)^(n>>j))&1)<<j))^n

log en base 2
log2(x)=log(x)/log(2)

decomposition en fractions continue du reel r en i étapes
e.g. cfraci(pi,3)= "355/113"
cfrac(r,i)=(r==floor(r))?"decimal svp":cfr3(r,i,floor(r),0,1,1,0)
cfr3(r,i,a,h2,k2,h1,k1)=(i<=0)?pfrac((a*h1+h2),(a*k1+k2)):cfr3(1./(r-a),i-1,floor(1./(r-a)),h1,k1,a*h1+h2,a*k1+k2)
pfrac(a,b)="".a."/".b

pareil que cfrac mais rend un complexe {a,b} au lieu de la chaine "a/b"
e.g. cfraci(pi,3)= {355.0, 113.0} c'est à dire 355+113i
cfraci(r,j)=(r==floor(r))?{0,0}:cfr4(r,j,floor(r),0,1,1,0)
cfr4(r,j,a,h2,k2,h1,k1)=(j<=0)?(a*h1+h2)+i*(a*k1+k2):cfr4(1./(r-a),j-1,floor(1./(r-a)),h1,k1,a*h1+h2,a*k1+k2)

min(a,b)=(a<b)?a:b
max(a,b)=(a>b)?a:b
frac(x) = x-floor(x)

affiche les valeurs de f(x) entre x=a et x=b inclus, par incrément de 1, sous forme d'une chaîne.
a et b doivent être entiers, avec a<=b
par exemple si f(x)=x**2-4 alors vals(1,5) = "-3,0,5,12,21,"
vals(a,b)=(b<a)?"":(b==a)?"".f(a):"".vals(a,(a+b)/2).",".vals((a+b)/2+1,b)

donne tous les entiers n entre a et b inclus tels que f(n)=0 .
par exemple si f(x)=x**2-4 alors zeros(-5,5) = "-2,2"
zeros(a,b)=(a>b)?"":(f(a)==0)?"".a.",".zeros(a+1,(a+b)/2).zeros((a+b)/2+1,b):zeros(a+1,(a+b)/2).zeros((a+b)/2+1,b)

pi est déjà prédéfinie. J'ai ajouté :
i = {0,1}
e = 2.71828182845905 # The base of the natural logarithm
sqrt2 = sqrt(2) # The square-root of 2
euler = 0.57721566490153 # The Euler-Mascheroni constant